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ABSTRACT 

Ring polymer molecular dynamics (RPMD) is used to directly simulate the dynamics of an excess electron in a 
supercritical fluid over a broad range of densities. The accuracy of the RPMD model is tested against numerically exact 
path integral statistics through the use of analytical continuation techniques. At low fluid densities, the RPMD model 
substantially underestimates the contribution of dclocalized states to the dynamics of the excess electron. However, 
with increasing solvent density, the RPMD model improves, nearly satisfying analytical continuation constraints at 
densities approaching those of typical liquids. In the high density regime, quantum dispersion substantially decreases 
the self-diffusion of the solvated electron. In this regime where the dynamics of the electron is strongly coupled to the 
dynamics of the atoms in the fluid, trajectories that can reveal diffusive motion of the electron are long in comparison 
to (3h. 
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I. INTRODUCTION 

Ring polymer molecular dynamics (RPMD) allows for the direct simulation of quantum mechanical systems. It is a 
classical model that both preserves the exact quantum Boltzmann distribution and exhibits time-reversal symmetryfi^ 
These features ensure that RPMD trajectories are stable and self-consistent on long timescales, enabling the study 
of coupled dynamical timescales in complex systems. In exhibiting these features, RPMD and centroid molecular 
dynamics 3 -" 4 - are unique among quantum dynamical methods. Mixed quantum-classical methods based on mean- 
field^ and trajectory surface hopping 7 ."^ dynamics do not preserve the correct Boltzmann populations 3 i 10 i n Similarly, 
methods based on the classical Wigner mode l 12 i 13 ' 14 i 15 ' 16 ' 17 i 18 ' 19 do not yield trajectories that preserve the quantum 
Boltzmann distribution and are thus not employed as models for direct dynamical simulation. 

Here, we test the RPMD model for the direct simulation of an excess electron in a supercritical fluid. This is the first 
application of RPMD to electronic degrees of freedom, a highly quantum mechanical - and presumably challenging - 
regime for the model. To evaluate the accuracy of RPMD in this context, we compare the model dynamics against 
numerically exact path integral statistics with the aid of analytical continuation techniques. We further explore the 
long-timescale features of the dynamics using direct simulation with RPMD. 



II. METHODS 

An excess electron is simulated in an otherwise classical fluid. The system is described using potentials and 
parameters that were previously adopted by Berne, Coker, and coworkers to model an excess electron in helium 
at 300 K i 20 i 21 i 22 i 23 i 24 i 25 i 26 ' 27 The potential energy function, U(q, Qi, . . . , Qat), is a sum of pairwise helium-helium 
and electron-helium interactions^ q and Qj specify the positions of the electron and the helium atoms, and N is 
the number of atoms in the fluid. Helium-helium interactions arc described using the Lennard- Jones potential with 
(7Hc = 2.556 A and ejjc = 10.22 K; electron-helium interactions are described using the pairwise pseudopotential 
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where r is the electron-helium distance and, in atomic units, A = 0.655, B ~ 89099, and C — 12608. 
The ring polymer molecular dynamics (RPMD) equations of motion for this system arei 
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where the integer a indexes the ring polymer beads for the electron, such that q(°) = q( n '. Also, V" and Vj specify 
the respective velocities for the ring polymer beads and the classical helium atoms, m and M are the respective 
electron and helium masses, and u) n = n(ph)^ 1 , where j3 is the reciprocal temperature. These dynamics preserve the 
path integral discretization for the Boltzmann distribution ; 28 ! 29 



P({q( Q )}, Qi, . , Q N ) oc exp j-£ (^n(q (a) ~ q (Q " 1} ) 2 + U(rf a \ Qi, Q*)) | 
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such that in the limit of large n, sampling the RPMD trajectories yields exact equilibrium averages . 30 ! 31 For example, 
the imaginary time mean-squared displacement 3 ^ 

i? 2 (r)H|q(T)-q(0)| 2 >, (4) 

is obtained from the RPMD using 

T 

R 2 (r a ) = ((qW - q<°)) 2 ) np = J- £ (qW(,) - ^(stf , (5) 
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where r Q = j/3h/n, and q( a )(t) indicates the position of bead a after evolution for time t using Eq. @. T b 8 is the 
total length of the RPMD trajectory, which must be thermostatted to be fully ergodic. For the calculation of static 
quantities, RPMD is equivalent to the path integral molecular dynamics (PIMD) method i 30 i 31 

Beyond calculating equilibrium properties, RPMD utilizes the time-displaced statistics to estimate real-time be- 
havior. It has been employed previously as a method to calculate Kubo-transformed correlations functions^ ' 33 ! 34 ' 35 ! 36 
such as the velocity autocorrelation function to which the mobility of the excess electron and its absorption spectrum 
can be related, 
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5v.v(t) = -30 / d\ (p(-*Aft) • p(t)>. (6) 
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Here, p = — ifiV q is the momentum operator for the electron. The RPMD approximation to this correlation function 

j : 33.34 

c v .v (t) « (v(0) • v(i)) Rp = J- g v( S ) • v(s + t), (7) 

iobS 8=1 

where v(t) indicates the bead-averaged momentum for the electron ring polymer that is also evolved according to Eq. 
©. 

In the current work, we take a somewhat expanded view of the RPMD model, using it to directly simulate the 
time-evolution of a quantum system. For example, rather than calculating the self-diffusion coefficient of the excess 
electron from the integral of c v . v (t), we instead simulate the real-time displacement of the electron using 

i? 2 (t)«((q(t)-q(0)) 2 ) Rp , (8) 



where q(t) is the bead-averaged position for the time-evolved electron ring polymer. The RPMD model obeys the 
classical molecular dynamics relationship 37 

lf t R 2 (t) = fdt'c v . v (t'), (9) 

so that both approaches yield the same estimate for the self-diffusion coefficient^ As another example, we note that 
the RPMD estimate for the thermal rate constant need not be obtained from the calculation of a time correlation 
function; an equivalent estimate would be obtained by simply running a long RPMD trajectory and counting the 
frequency with which the reaction of interest occurs! 3 -* 

At the heart of these and other robust features of RPMD is the status of the model as a genuine classical dynamics, 
albeit one that preserves the quantum Boltzmann distribution. RPMD ensures time reversibility, and it exhibits time- 
translational invariance at equilibrium. It is a model for simulating the real-time evolution of a system, in addition 
to a means of calculating time correlation functions. RPMD can thus be employed outside of the time correlation 
function formalism and beyond the linear response regime. In the same way that standard molecular dynamics is 
used to simulate classical reaction mechanisms and to simulate classical processes far from equilibrium, we propose 
that RPMD be used to simulate dynamics in systems for which quantum statistical effects play an important role. 
The centroid molecular dynamics model shares these features and prospects, although it has not, to our knowledge, 
fully utilized this philosophy. 

A. Relating static and real-time correlation functions via analytical continuation 

To evaluate the accuracy of RPMD for the dynamics of the solvated electron, we will relate the results of the model 
dynamics to numerically exact path integral statistics through the use of analytical continuation ! i This approach 
is described in the next two sections. 

The dipole autocorrelation function for the solvated electron is 

C(t) = (n(t) • MO)} = ^Ti{e- f>& e it& / h iuT itk l n ■ n). (10) 

Here, Z = Ti(e~^ H ) is the canonical partition function, fi — — |e|q is the dipole operator, e is the charge of an electron, 
and (3 = 1/(&bT) is the reciprocal temperature. The dipole spectral density function, I(uj), is defined such that 

/OO -t POO 

dt e lult C{t) and C(t) = — / du e~ lut I{uj). (11) 
.~, 27T . _^ 



The correlation function C(t — it) is analytic in the strip < r < f3h, where t and r are real numbers. We can thus 
introduce the imaginary-time correlation function G(r) = C(— it) such that, by way of analytic continuation, G(r) 
formally contains all of the dynamical information on the real-time axis4^ 

It follows from Eq. (fTTj) that G(r) and I(oj) are related by a double-sided Laplace transform, 

G(t) = I duj e- UT I(uj). (12) 

By discretizing G(r) in imaginary time, it is then straightforward to show that 

R 2 (r) = |(G(0)-G(r)) 

= -L I" du>I(u>)(l-e- u T), (13) 

where R 2 (t) is the mean squared displacement of the electron in imaginary time from Eq. Finally, using the 

detailed balance relation I(—u>) — e~ P huJ I [lu) , we obtain 

fir r°° 

R2 (r) = ink / fa k(t,u)*{l,), (14) 
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where c is the speed of light. Here, we have introduced the absorption cross section 
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Eqs. (fbI|) - (flT)|) were previously reported by Gallicchio and Bcrnc 
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B. Maximum Entropy Analytical Continuation 

The analytical continuation result in Eq. (|14[) suggests an integral inversion problem that yields real-time dynamics 
(a(to)) from purely statistical information (R 2 (t)). However, this inversion is known to be numerically unstable, 
such that small errors in R 2 (t) lead to large errors in cr(uj). Nonetheless, the inversion can be performed indirectly 
with the aid of information theory methods such as the maximum entropy technique i 41 i 42 ' 43 ' 44 Maximum entropy 
analytical continuation (MEAC) yields cr(w) by refining a prior estimate for the spectrum, a°(uj), against a numerical 
calculation for R 2 (t). In the current study, we shall obtain <j°{lo) from the RPMD model, and we shall obtain R 2 {t) 
from numerically exact path integral statistics. The degree to which the maximum entropy inversion alters the RPMD 



prior provides insight to the accuracy of the RPMD model for the excess electron. The idea of using MEAC to correct 
approximate RPMD correlation functions was recently put forward by Habershon et al. , and MEAC was previously 
applied to combine centroid molecular dynamics with imaginary time data4£ Our implementation of the maximum 
entropy inversion follows that of Habershon et alM^ and Gallicchio and Berne ; 26 ! 27 which are in turn based on the 
prescription provided by Bryan 4£ 

In our numerical implementation, R 2 (t) is known on the discrete set of n points {tj}, and we wish to solve for a(ui) 
on the discrete set of M points {oJk}- Eq. (fT4|) is thus expressed as a matrix equation 

c = Ks (17) 

where c is a column vector with elements Cj = R 2 (jj), K is a rectangular matrix with elements Kjk = K(jj, u>k), s is 
a column matrix with elements Sk — cr(wfe)Awfe, and Aw^ is a quadrature weight. The maximum entropy approach 
does not seek to directly invert this equation, but rather to maximize the scoring function 

Q(s;a) = aS(s)- X 2 (s)/2 (18) 

with respect to the elements of s. In this equation, % 2 (s) is a measure of the degree to which s fits the imaginary time 
data, S(s) is a measure of the degree to which s deviates from an initial guess for the solution, and a is a parameter 
that tunes the relative influence of these measures during maximization. 
Specifically, 

X 2 (s) = (c-Ks) T C- 1 (c-Ks), (19) 

where C is the covariance matrix 

1 M 

^-Aic^£( w -^)(fe>-^)- (20) 

and M. is the number of statistically independent measurements of the imaginary time correlation function R 2 (jj) 
obtained during the equilibrium sampling of the quantum Boltzmann distribution. 
Given a prior model for the spectrum cr°(a->fc), we can define the information entropy 
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where — a (uk)AoJk- To obtain the RPMD prior spectrum a°{uj), we note that absorption cross section satisfies 

, , 47re 2 /3 ; , , 
<r(w) = -Cv. v (w), (22) 



where c v . v (ui) is the Fourier transform of the Kubo-transformed velocity autocorrelation function. The RPMD model 
for the time correlation in Eq. ([7|) thus yields 

a°(u) = £ / dt e~ luJt c v . v (t). (23) 



Finally, the parameter a in Eq. (fl8|) was obtained using the L-curve method.^'- - In this method, Q(s;a) is 
maximized for a series of a. For each value of a, the quantity ln[% 2 (s)] is plotted against ln[— 5(s)], and the resulting 
curve exhibits a characteristic "L" shape. The bend in this curve marks the transition from a regime in which 
deviations from the prior model are rewarded with substantial improvements in the degree to which the solution fits 
the imaginary time data, to a regime in which further deviations from the prior model yield smaller improvements in 
the fit of the imaginary time data. Results are reported at the value of a that coincides with this transition (i.e., the 
point of maximum curvature in the L-curve). 

We emphasize that analytical continuation can only provide a quality check for the short-time RPMD dynamics, 
because it offers little information about timescales that exceed [3h. This well known limitation of analytical continu- 
ation can be illustrated by considering the degree to which different dynamical timescales are weighted by the kernel 
K (u>, t) in Eq. (fT4"|) . For < t < f3h and t > 0, the Fourier transform of the kernel is 

/oo 
dt e- iuJt K(u,T) 
-OO 

This result, which is obtained by straightforward contour integration after observing that the K(lu,t) exhibits a string 
of simple poles for purely imaginary values of to, shows that the degree to which the kernel reports on large t vanishes 
exponentially. For each value of t, the largest contribution to K(t,r) occurs at r = 0H/2, where we obtain the closed 
form expression 



K(t,0h/2) = 2 In 



coth | ) 



(25) 



This expression decays exponentially at long times with a rate of (3h/(2i:). The same rate for exponential decay was 
reported by Habershon et al. for analytical continuation involving a different kernel^ 



C. Simulation details 



All simulations are performed using N — 1000 classical helium atoms in a cubic simulation cell with periodic 
boundary conditions. The electron is represented using n = 1024 ring polymer beads. The simulation temperature is 
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309 K, and the helium fluid is studied at reduced densities of p* =0.1, 0.3, 0.5, 0.7, and 0.9. The values of n and N that 
we employ are slightly larger, but comparable, to the values used in previous path integral Monte Carlo simulations 
of this system under the same conditions. 22 ! 26 ' 27 The classical limit for the dynamics of the electron is obtained by 
setting n = 1. 

RPMD trajectories are performed by integrating Eq. ([2]). Because of the range of timescales in the problem, 
some care is needed to do this integration efficiently. The multiple time-stepping (MTS) algorithm is employed to 
overcome the large difference in timescale for the electronic and helium degrees of freedom i 50 ' 51 In all simulations, the 
classical helium atoms are evolved with a timestep of 0.32 fs; during each of these nuclear timesteps, the coordinates 
and forces of the ring polymer for the electron are updated 170 times. As in previous RPMD simulations, each 
timestep for the electron ring polymer involves coordinate updates due to (1) forces from the physical potential 
(— Vq(*o U (q^ k \ Qi, • • ■ , Qjv)) and (2) the exact evolution of the purely harmonic portion of the ring polymer potential. 
The combined integration scheme is time-reversible and symplectic. 52 

All interactions are truncated at a cutoff distance of 2.5 ctho- Because of our MTS integration scheme, the primary 
computational expense arises from electron-helium force evaluations. However, these force evaluations are easily 
parallelized to yield a considerable (~20-fold) speedup in the wall-clock time needed for the simulations. 

Prior to performing the recorded RPMD trajectories, the system is equilibrated at each density as follows. First, 
the helium fluid is equilibrated in the absence of the excess electron using a long, thermostated molecular dynamics 
trajectory at 309 K. Then, the electron ring polymer is inserted into the equilibrated configuration of the helium 
fluid, and the combined system is run for 200 RPMD trajectories of length 380 fs. Between trajectories, the position 
coordinates for the system are not changed, but the momenta are resampled from the Maxwell-Boltzmann distribution. 
These equilibration RPMD trajectories are discarded. 

The recorded RPMD trajectories are continued from the end of the equilibration trajectories. At each density, we 
perform 30 recorded trajectories of length 7.66 ps. Each recorded trajectory includes nt = 24000 MTS integration 
timesteps. As before, the momentum coordinates are resampled between trajectories. Eqs. (|20|) . ([7]), ([8]), and ([5]) are 
employed during these recorded trajectories to sample the quantities of interest. 

The RPMD results for R 2 (t), Cij, and u°{lo) provide the input to the MEAC calculation. The imaginary time 
correlation function R 2 (t) is known on the discrete set of points {tj}, where Tj — ^f3h and j = 1, 2, . . . , n. Likewise, 
the covariance matrix Cy for the imaginary time correlation function is obtained on i, j = 1, 2, . . . , n. Using Eq. (|23[) . 
the RPMD prior spectrum cr°(w) is obtained by Fourier transforming the RPMD velocity autocorrelation function 



9 

c v . v (t) on the discrete set of points {wfc}, where u>k = ^w max , S.w max = 6.25 eV, k = 0, 1, . . . ,7V, and J\f — n t /2. 

Unphysical noise can be eliminated from a spectrum by damping the corresponding time correlation function at 
large times before performing the discrete Fourier transform. Thus, in preparing <t°(w) for the MEAC calculation, we 
damp c v . v (i) at times greater than id using the weighting function 

w(t) =cxp[-(t-t d ) 2 /(2((3h)% (26) 

where t d = 2j3h « 50 fs for p* = 0.1, and t d = (3H w 25 fs for p* = 0.3 - 0.9. The MEAC calculations were also 
performed using larger values of t d to confirm that our conclusions are independent of this parameter. 

III. RESULTS AND DISCUSSION 
A. Dynamics on short timescales 

Fig. [T]A_ shows the imaginary time correlation function R 2 (t) that is sampled during the RPMD trajectories at 
various densities using Eq. ([5]). These numerically exact results are equivalent to previously reported path integral 
Monte Carlo simulations for the system considered here , 22 ' 26 ' 27 and they illustrate the finite-temperature manifestation 
of Anderson localization for an excess electron in a disordered mediunii 53 i 54 

For a given solvent configuration, the electron energy eigenfunctions are either spatially extended, in which they 
percolate across the periodic simulation cell, or they are spatially localized within the celli 2 ^ The relative stability of 
extended vs. localized wavefunctions hinges on a quantum mechanical balance between the penalty for delocalizing 
the electron wavefunction in a quasi-disordered medium and the penalty for creating a solvent cavity that is large 
enough to confine the localized electron.— This balance shifts as a function of fluid density. 

Since R 2 ((3h/2) is a measure of the spatial extent of the excess electron, the trend of increased localization with 
increasing fluid density is clearly seen in Fig. [T]A For p* = 0.1, the environment of the electron is only weakly perturbed 
from that of a free particle. The fluid assumes configurations in which thermally accessible electron eigenstates are 
delocalized, giving rise to an imaginary time correlation function that approaches the free particle result, 

Rl cc (r)=3\ 2 T(ph-T)/(ph) 2 , (27) 

where A = (fi^fi/m) 1 / 2 is the thermal wavelength of the electron. As p* increases, the atoms in the fluid create a more 
densely disordered environment that destabilizes spatial coherence in the electron eigenfunctions. The fraction of 
thermally accessible solvent configurations with thermally accessible delocalized electronic states diminishes, and the 
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system experiences a thermodynamic driving force to localize the electron in a solvent cavity, i.e. to form a polaron. 

Fig. [TJ3 shows the velocity autocorrelation functions c v . v (t) that are obtained from direct RPMD simulations of 
the excess electron using Eq. ([7]). Like the imaginary time correlation functions, these real-time correlation functions 
depend strongly on the solvent density. For p* = 0.1, c v . v (t) decays relatively slowly and exhibits only a weak rebound 
due to collisions with the dilute solvent environment. However, with increasing helium density, the ballistic motion 
of the electron in the RPMD model decays much more rapidly. The timescale for the first collisional rebound of the 
electron (i.e., the first minima of c v . v (t)) decreases from 7 fs at p* = 0.1 to 1 fs at p* — 0.9. Furthermore, as the 
density approaches p* = 0.9, we see the emergence of multiple peaks in the correlation function, or "rattling " of the 
excess electron in its solvent cage. 

In Fig. [Tp, we report the RPMD estimates for the absorption cross section of the solvated electron <t°(lu), which 
are obtained from the c v . v (t) using Eq. (|23[) . With increasing helium density, the absorption peak is shifted to 
higher frequency, indicating an increasing energy gap between the ground-state and lowest-energy excited states 
of the electron. Greater solvent density also broadens the spectra and elongates the high-frequency spectral tail. 
All of these features are consistent with previous reports using the RISM-polaron theory for the excess electron in 
an adiabatic hard-sphere fluid 5 - and using numerical analytical continuation calculations for the system considered 
herei 26 ' 27 The only qualitative difference between the RPMD absorption cross sections in Fig. QTJ and those of previous 
theoretical studie a 27 ' 55 occurs in the low density regime, at p* =0.1. The RPMD model predicts a maximum in the 
absorption cross section at finite frequency, whereas previous calculations have found that the low-density spectrum 
decays monotonically from a maximum at u> = Q. 

Eq. (|14p shows that the imaginary time correlation function places a constraint on the dynamics of the solvated 
electron. The accuracy of the RPMD model is thus indicated by the degree to which it violates this dynamical 
constraint. In Fig. [TjD, we plot the RPMD prior estimate for the imaginary time correlation function, which is 
obtained by inserting a°(ui) into the right-hand side of Eq. (|14p . If RPMD provided an exact description of the 
solvated electron dynamics, then Fig. [T}D would precisely mirror the exact imaginary time correlation functions in 
Fig. [T]A_. In the high density regime, this is very nearly the case; the RPMD model dynamics are consistent with the 
imaginary time data for p* > 0.5. However, at p* = 0.1, and to a lesser extent at p* = 0.3, the RPMD estimate 
deviates substantially from the exact result. In this low-density regime, the RPMD prior spectrum (J°{lo) corresponds 
to the dynamics of a model electron that is more localized than the actual solvated electron. 
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Since the imaginary time correlation function can be written in the basis of electronic eigenstates as follows, 



we see that the changes in the value of this function with t arise from the accessibility of excited electronic states. 
At low helium densities, the fact that the imaginary time correlation function for the RPMD prior model (Fig. [Tp) 
plateaus at smaller values of (r/(3H)(l — t/ ' fiK) than the numerically exact result (Fig. [T]A_) suggests that the RPMD 
model underestimates the role of these excited states in the dynamics of the electron. By the same argument, the 
dominance of the electronic ground state in the dynamics at higher helium densities is correctly captured by the 



Eq. (|Ti)) can also be used to obtain a first-order correction to the RPMD prior spectrum u°(uj). We employ the 
maximum entropy analytical continuation (MEAC) technique to numerically invert Eq. (| 14[) . using (t°(oj) as an initial 
guess for the solution cj(lo). As was explained in Sec. Ill Bl this approach alters the RPMD prior spectrum to yield a 
solution that reproduces the imaginary time correlation functions in Fig. [TJ\ to within statistical certainty. 

Fig. [U presents the MEAC correction to the RPMD model at each density. The dashed lines in the left column 
show the RPMD correlation functions c v . v (t), and the dashed lines in the right column show the RPMD prior spectra 
<7°(w). These results are unchanged from Figs. [TJ3 and[Tp, respectively. The solid lines in the right column present 
the MEAC corrections to the RPMD prior spectra. Finally, the solid lines in the left column show the velocity 
autocorrelation functions that correspond to the MEAC-corrected spectra. 

Both columns of Fig. [2] illustrate that the MEAC procedure alters the RPMD model less with increasing helium 
density. This correction generally narrows, red-shifts, and intensifies the absorption band for the RPMD prior. For 
p* = 0.1, the MEAC correction is most dramatic, causing a substantial increase in the absorption intensity and a 
shift in the absorption peak to nearly ui = 0. This result again suggests that the RPMD model underestimates the 
role of low-lying electronic excited states in the dynamics of the electron at low solvent densities. For p* > 0.3, the 
MEAC correction is much reduced, and it no longer changes the qualitative features of the RPMD prior. Indeed, for 
p* = 0.7 and 0.9, the RPMD time correlation function and the MEAC-corrected results are very similar. 

It is not surprising that the RPMD model is more reliable at high densities than at low densities. In the low density 
regime, where the electron achieves high mobility via thermal access to extended electronic states, the dynamics of 
the electron is a coherent scattering problem that is beyond the realm of applicability for the simple RPMD model. 
On the other hand, at higher densities, the dynamics of the electron is almost entirely ground-state dominated and 
governed by the electronically adiabatic diffusion of the solvent-electron polaron. In this case, since RPMD exactly 




(28) 



RPMD model. 
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describes both the structural (i.e. static) features of the localized electron and the classical dynamics of the helium 
solvent, it is expected that the RPMD model is accurate. From Fig. [21 we are encouraged that the MEAC correction 
to the RPMD model is minor for the densities of typical liquids. 

The MEAC-corrected spectra reported in the right column of Fig. [5] can be compared with previous MEAC studies 
for the same system.! 26 ' 27 The primary way in which our calculations differ from this earlier work is in the choice of 
prior spectrum for the maximum entropy inversion' we use the RPMD model to obtain the prior spectrum, whereas 



Rcf. 



26l employed a flat prior spectrum and Ref . 



271 employed a prior spectrum that is based on sum rule constraints 
on the imaginary time correlation function. Differences in the results of these three studies illustrate that, given the 
finite statistical error in the calculation of the imaginary time correlation function, the maximum entropy inversion 
of Eq. rfj] is unable to fully overcome different choices for the prior spectrum. Nonetheless, we find that the MEAC 
technique is a useful means of testing the short timescale description of simple dynamical models such as RPMD. 



B. Dynamics on long timescales 

The left column of Fig. [3] presents the mean squared displacement (MSD) for the electron using the RPMD model 
(dashed lines) for times up to 200 fs. For comparison, we also include the MSD for the atoms in the helium solvent 
(solid line), and at higher densities, we plot the linear fit to the RPMD curve at t = 25 fs, (i.e., t = f3h). It is clear, 
especially at higher densities, that the dynamics of the electron does not reach the diffusive regime on the timescale of 
/3H. The MSD curves continue to deviate from linear behavior for many times this thermal timescale. This is at first 
surprising, since the RPMD velocity autocorrelation functions appear to decay quickly in Fig. [21 but the apparent 
contradiction simply illustrates that the self-diffusion of the electron is coupled to the slower dynamics of the solvent 
atoms. This is easily understood in the high density limit, where the delocalized states of the electron are thermally 
inaccessible from the strongly localized ground state and the mobility of the electron arises almost exclusively from 
the non-adiabatic diffusion of the solvent-electron polaron. In this limit, the mobility of the electron would approach 
exactly zero if the solvent were fixed^ so any observed diffusion of the electron must be coupled to movements in 
the solvent environment. For all p* > 0.5, the left column of Fig. [3] suggests that coupling of the electronic and 
solvent dynamics plays an important role. The rapid decay of the velocity autocorrelation functions in Fig. [2] is thus 
a graphical illusion; as is required by the mathematical equality in Eq. @, integration of the long-time tail reveals 
its non-zero contribution to the slope of the mean-squared displacement. 

The importance of dynamical timescales that are slow in comparison to (3H also deserves commentary from a 
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methodological perspective. As we have already pointed out in connection to Eq. ([25]) , the MEAC technique contains 
exponentially little information about dynamics on these long timescales, which explains why previous theoretical 
studies of the solvated electron have not observed the slow onset of diffusive behavio n 26 ' 27 Furthermore, the impor- 
tance of long timescale dynamics underscores the need for a model dynamics that rigorously preserves the quantum 
Boltzmann distribution, such as RPMD or centroid molecular dynamics. For methods such as the linearized semi- 
classical initial value representation (SCTVR) ) 12 ' 13 ' 19 in which classical molecular dynamics trajectories are initialized 
from a quantum distribution, it is expected (because of the negligible mass of the electron in comparison to the com- 
bined mass of the solvent) that the long-time asymptote of the electron dynamics will be independent of the initial 
distribution and will incorrectly converge to a purely classical description. Conclusions about the coupled dynamics 
of the electron and solvent may then be obscured by the gradual breakdown of the statistical description with time^ 
The right column of Fig. [3] presents the RPMD results for the electron MSD on the timescale of picoseconds (i.e., 
up to w 150/3/i). Again, we include the solvent MSD for comparison. The straight lines in the left column of Fig. [3] 
indicate that both the electron and the solvent have reached the diffusive regime. The slope of the dashed line can 
be used to obtain the self diffusion coefficient for the RPMD model of the electron dynamics, 

D RP = i lim 4R 2 (t). (29) 
6 t-+oo at 

Eq. 12]) ensures that this expression for -Drp is equivalent to evaluating the integral of the RPMD velocity autocor- 
relation functions. 

In Table HI we report the RPMD electron self-diffusion coeffient Drp and the helium solvent diffusion coefficients 
Z?Hc- Also included is the self-diffusion coefficient D c \ for the electron obtain from classical molecular dynamics 
simulations (i.e., 1-bead RPMD). 57 Each is obtained from a linear least-squares fit of the corresponding MSD curve 
between 1 ps and 3.7 ps. Although both the classica l 58 ' 59 and the RPMD 33 ' 34 self-diffusion coefficients are subject 
to a hydrodynamic system size effect, it is not likely to alter our conclusions here^ Also in Table Q] we include an 
estimate for the radius of the excess electron at the various solvent densities. The size of the repulsive core for the 
electron is obtained by evaluating the distance at which the electron-helium interaction potential in Eq. |T]) first 
goes to zero, namely i?o = (B — C) 1 / 6 w 1.35<7H e - The effective radius of the classical model for the electron is thus 
Rci = Rq — 0.5<7hc ~ 0.85(7hc- In the RPMD model for the electron, the quantum dispersion of the electron swells its 
size by an amount that can be estimated from the radius of gyration of the ring polymer, such that Rrp = R c \ + R s , 
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where 




RP 



(30) 



Comparison of -Drp and D c \ in Table U reveals that at every density considered, inclusion of quantum effects with 
the RPMD model substantially decreases the mobility of the excess electron. This differs from the quantum effect 
on self-diffusion in more weakly quantum mechanical systems, such as bulk liquid wateri^i In a weakly quantum 
mechanical system, quantum dispersion softens the interaction potentials without greatly affecting the structure of 
the system or the mechanism of diffusion; the result is an enhanced rate of self-diffusion. The highly quantized excess 
electron, however, presents a much larger collisional cross section to the fluid atoms than does the corresponding 
classical model, resulting in greater drag and a lower self-diffusion coefficient. 

It is also interesting to compare the self-diffusion coefficients for the electron and the solvent atoms in Table U At 
the lowest fluid density, the RPMD model electron diffuses much more rapidly than the solvent atoms, but the reverse 
is found at all higher densities. To understand these results, we can compare the ratio of diffusion coefficients from 
our simulated dynamics, -Drp / -Due, with the Stokes-Einstein prediction for this ratio, 



where -Rrp is the calculated size of the quantized electron in Table [J and R^ c = 0.5crHe is the van der Waals radius 
of the solvent atoms. Eq. (|3ip simplistically assumes that the electron undergoes overdamped diffusion as a spherical 
solute of fixed radius -Rrp. 

In Fig. 21 we plot Drp / -Drc from the simulation results in Table U (solid) and from the Stokes-Einstein model 
in Eq. (f3"Tj) (dashed). At high densities, the different estimates agree, which is reasonable given that the electron 
is strongly localized in this regime. However, for lower densities, the diffusion of the electron in the RPMD model 
is enhanced relative to the prediction of the adiabatic Stokes-Einstein model. This deviation from Stokes-Einstein 
behavior is also reasonable because of the important role of non-adiabatic transitions to delocalized electronic states 
in the diffusion of the excess electron at low density. Although it was found in connection with Fig. Q] that RPMD 
underestimates the contribution of delocalized states to the mobility of the excess electron at low densities, we see 
here that the RPMD model at least partially captures this important effect. 




(31) 
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IV. CONCLUSIONS 

We have used ring polymer molecular dynamics to simulate the diffusion of an excess electron in supercritical 
helium at various densities. The RPMD model uses imaginary time path integrals to represent both the electron and 
the solvent atoms in the position basis, and it prescribes a classical molecular dynamics that rigorously samples the 
quantum Boltzmann distribution. In addition to conveniently placing the electronic and solvent degrees of freedom 
on the same dynamical footing, the RPMD model can be used to run long trajectories and to study the coupling of 
different dynamical timescales. 

With the aide of analytical continuation relationships, we have tested the dynamics of the RPMD model against 
numerically exact static correlation functions (Figs. [TJA. and [If), Fig. [2|. These tests, which are most stringent for 
timescales shorter than (3h, suggest that RPMD is increasingly accurate as the solvent density increases. Indeed, for 
p* > 0.5, the dynamics convincingly satisfies the analytical continuation constraint. This result is encouraging for the 
prospect of using RPMD to study an excess electron in liquid water. A simple estimate based on the peaks in the 
oxygen-oxygen radial distribution for ambient water reveals that that the liquid has an effective reduced density of 
p* « 0.5 — 0.9, a regime in which we find the RPMD model to be reliable. 

At lower solvent densities, our analytical continuation results indicate that the RPMD model systematically under- 
estimates the degree to which delocalized elecronic states enhance the mobility of the excess electron. However, the 
failure at low densities is not absolute. By comparing the -D^p/Due ratio from simulation with an estimate based 
on the Stokes-Einstein diffusion of the excess electron, we find that the RPMD model does predict a substantial 
enhancement in the electron diffusion at low densities. Even though RPMD does not fully account for the increase in 
electron mobility due to thermally accessible delocalized electronic states at low densities, it does qualitatively include 
the effect. 

Finally, by running RPMD trajectories on the timescale of picoseconds, we have been able to fully bridge the 
timescales of the electronic and the solvent atom dynamics. In doing so, we find that these timescales become coupled 
as the electron localizes with increased solvent density. This fundamental example of complex dynamics is conveniently 
probed using the RPMD model, and it is clear that RPMD can be similarly employed to study the coupled dynamics 
of quantum mechanical and classical mechanical motions in other contexts. 
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TABLES 



TABLE I: Self-diffusion coefficients" for the electron and the solvent helium atoms and the radius of the solvated electron 6 at 
various densities. 



p* 


I>RP 


Dei 




Rbp 


0.1 


1.1(1) 


6.6(3) 


0.207(4) 


3.69 


0.3 


0.032(3) 


0.53(3) 


0.0702(6) 


2.71 


0.5 


0.014(1) 


0.056(3) 


0.0391(4) 


2.26 


0.7 


0.0085(9) 


0.022(1) 


0.0256(4) 


2.03 


0.9 


0.0043(4) 


0.0119(6) 


0.0176(2) 


1.87 



"All self-diffusion coefficients are reported in units of ag/fs. The uncertainty in the final digit is indicated in parenthesis. 
b The radius of the quantized electron, i?R,p, is reported in <7 He - The radius of the electron in the classical model is 
R c i = 0.85(7hc, and the radius of the solvent helium atoms is 0.5<th c - See text for details. 



20 



FIGURE CAPTIONS 



Figure 1. Data from the RPMD simulations of the excess electron in helium at various densities. (A) Numerically 
exact imaginary time correlation functions R 2 {t). (B) RPMD velocity autocorrelation functions c v . v (t). (C) RPMD 
prior spectra <j°{u>). (D) RPMD prior estimates for the imaginary time correlation functions, which are obtained by 
inserting the cr°(w) into Eq. (|14[) . The imaginary time correlation functions in parts A and D are symmetric about 

t = fih/2. 

Figure 2. The ME AC correction to the RPMD model at various densities. The left column presents the RPMD 
(dashed) and MEAC-corrected (solid) velocity autocorrelation functions for the electron. The right column presents 
the corresponding absorption cross sections. 

Figure 3. The RPMD mean squared displacments for the electron (dashed) and the solvent helium atoms (solid) at 
various densities. Results are presented for both short timescales (left column) and long timescales (right column). 
Figure 4. Ratio of the electron and helium atom self-diffusion coefficients at various densities, obtained from the 
RPMD simulations (solid) and from the Stokes-Einstein model in Eq. (|3Tj) (dashed). 
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